Hasta el momento hemos trabajado independientemente con caminantes aleatorios y distribuciones de probabilidad. Ahora trataremos de juntar los dos tópicos.
Si $P_{t}(i)$ es la probabilidad de que un caminante se encuentra en el sitio $i$ al tiempo $t$, entonces la distribución de probabilidad está dada por el conjunto $\{ P_{t}(i): i \in \mathbb{Z} \}$. Abstrayendo un poco, este objeto puede ser visto también como un vector con el número de entradas igual al número de sitios en el cual puede estar nuestro caminante.
A este vector lo llamaremos $\mathbf{P}_{t}$.
Nota 1: Este es el primer ejemplo en cual podemos decir que el tiempo es discreto, por lo que $t \in \mathbb{N}$.
Nota 2: En principio, el lugar donde nuestros caminantes marchen puede ser infinito (por ahora en una dimensión), por lo que el vector $\mathbf{P}_{t}$ tendría una infinidad de entradas.
Coloquemonos en nuestro paso $0$, i.e. $t=0$. El caminante estará en su condición inicial esperando a la bandera de salida. Supongamos que su condición inicial es en $i = 0$ y que el caminante tiene una probabilidad $p$ de dar un paso a la derecha y $q:= 1-p$ de darlo a la izquierda.
Para lo siguiente supondremos que $p = q = \frac{1}{2}$
No es dificil llegar entonces a que $P_{0}(i) = 0, \ \forall i \neq 0$ y que $P_{0}(0) = 1 $
En el primer paso, tendremos que: $$ \begin{matrix} P_{1}(-1) = \frac{1}{2} & P_{1}(0) = 0 & P_{1}(1) = \frac{1}{2} \end{matrix} $$
¿Qué pasa en el siguiente paso?
El espacio que el caminante puede abarcar se hace más grande, yendo de $i =-2,\dotsc, 2$. Ahora veamos como queda la distribución de probabilidad.
$$ \begin{matrix} P_{2}(-2) = \frac{1}{4} & P_{2}(-1) = 0 & P_{2}(0) = \frac{1}{2} & P_{2}(1) = 0 & P_{2}(2) = \frac{1}{4} \end{matrix} $$Las cosas ya se han puesto un poco más interesantes. Para entender un poco más como se han calculado un poco más rigurosamente estas probabilidades tomemos el caso de $P_{2}(0)$.
en el paso $t = 1$, el caminante tenía $\frac{1}{2}$ de probabilidad de estar en $i= -1, 1$. Supongamos que estaba en la celda $i=-1$, entonces el caminante tiene $\frac{1}{2}$ de probabilidad de dar un paso a la derecha en el paso $t = 2$; y de la misma forma, si el caminante estuviera en la celda $i = 1$, también habría $\frac{1}{2}$ de probabilidad de que al siguiente paso estuviera de nuevo en la celda $i = 0$.
De esta manera llegamos a la ecuación maestra de nuestro ejemplo en una dimensión y con probabilidades $p =q = \frac{1}{2}$
$$P_{t+1}(i) = \frac{1}{2}P_t(i-1) + \frac{1}{2}P_t(i+1) $$La generalización de la ecuación maestra es la siguiente
$$ \begin{equation} P_{t+1}(i) = pP_t(i-1) + (1-p)P_t(i+1) \ \ \ \ \ (1) \end{equation} $$[1a] Pongamos manos a la obra. El objetivo general será graficar cómo evoluciona temporalmente la distribución de probabilidad con ayuda de nuestra ecuación maestra.
Para esto necesitaremos un kernel un poco más sofisticado de los que ya hemos hecho. El gran cambio es que aquí sólo jugaremos con un solo caminante aleatorio, y no con varios como lo habíamos hecho hasta entonces. La malla de bloques de núcleos pasa a ser el espacio en que se mueve el caminante. Cambio sutil pero de grandes consecuencias.
En primer lugar necesitamos el arreglo en el cual el caminante pueda moverse de un lado a otro. A este le llamaremos $X$. Recuerda hacerlo lo suficientemente grande para que el caminante no choque con los extremos.
Para calcular la distribución de probabilidades del paso $t+1$ necesitamos la distribución del paso $t$. Sin embargo al sobreescribir nuestro arreglo $X$ perderemos información, y por lo tanto nuestros cálculos serán incorrectos. Es por esto que necesitaremos declarar otro arreglo en el cual podamos copiar nuestra información del tiempo $t$ para calcular la distribución deseada.
Ahora, aquí es donde viene lo interesante. La manera de copiar los datos. Para esto nos basaremos en el ejemplo de tiled programmation hecho en el Notebook 6 de la primera parte de este curso el cual se basaba en declarar un arreglo tipo __shared__ en el cual copiar los datos.
Así que veamos un poco como se vería el kernel. Supongamos que, tal cual dicho anteriormente, nuestros datos estuviesen en un arreglo $X$ y el lugar en el cual nos apoyaremos sea un arreglo en la memoria compartida llamado $X_{copia}$.
La idea general es entonces es de ir calculando los estados de X en el tiempo $t+1$ tesela por tesela. Supongamos entonces que X consiste de un arreglo de 200 celdas y buscamos calcular estas 200 celdas en el tiempo $t+1$ en grupos de 5. Entonces en $X_{copia}$ habremos de tener todos aquellas celdas con las cuales podamos realizar dichos cálculos.
En este caso en específico, puesto que el estado de una celda en el tiempo $t+1$ está determinada por ella misma y sus vecinos, entonces habremos de copiar cada uno de estos en $X_{copia}$. Esto no causa problemas a no ser por las celdas extremas de cada bloque. Para resolver esto tendremos que copiar también los vecinos que no aparecen en nuestro bloque de 5 celdas pero que también son necesarios para los cálculos.
De esta manera, para un arreglo de 200 celdas cuyos estados quieren ser calculados en bloques de 5, entonces necesitaremos en la memoria compartida teselas de 7 celdas.
A continuación mostramos la manera en la cual se copian los datos a la memoria compartida. En primer lugar mostraremos un programa en Python para darnos una idea más clara de que es lo que estamos buscando. Sólo después pasaremos al kernel.
In [2]:
import numpy as np
Supongamos un arreglo A con 17 estados iniciales. Nuestra intención es calcular el estado de cada celda en el tiempo siguiente. Usaremos el método con teselas para copiar los datos. Cada tesela se ocupará de 4 datos en A, por lo que según la ecuación maestra (1) necesitaremos que la dimensión teselar sea de 6.
In [3]:
A = np.array([1,2,3,4,5,6,7,8,9,10,11,12, 13, 14, 15, 16, 17])
tesela_A = np.ones(6)
In [9]:
# blockDim es la número de celdas que serán copiadas a la tesela
# gridDim es el número de bloques que tendremos
# ANCHO_TESELA es el número de datos que necesitaremos para calcular blockDim celdas en el siguiente estado
blockDim = 4
gridDim = len(A)/blockDim+1
ANCHO_TESELA = blockDim+2
# Regresamos a tener dos bucles...
# El primer bucle va de bloque en bloque por A
for blockIdx in xrange(gridDim):
# el segundo bucle busca los elementos de cada bloque
for tx in xrange(ANCHO_TESELA-1):
# los elementos necesarios son copiados a la tesela
if blockDim*blockIdx + tx-1 >= 0 and blockDim*blockIdx + tx-1 < len(A) :
tesela_A[tx] = A[blockDim*blockIdx + tx-1]
# y si hemos llegado a los extremos de A, entonces se coloca un 0.
else:
tesela_A[tx] = 0.0
# Este if else se dedica a colocar aquellos datos en la frontera derecha de la tesela
if 4*(blockIdx+1) < len(A):
tesela_A[ANCHO_TESELA-1] = A[blockDim*(blockIdx+1)]
else:
tesela_A[ANCHO_TESELA-1] = 0.
print B
Nota como en cada tesela los 4 valores centrales corresponden a aquellas celdas cuyo estado será calculado al tiempo $t+1$. En el caso de la última tesela, puesto que los valores de A ya fueron cubiertos, entonces la tesela es llenada con $0$'s para que no haya cálculos erróneos.
Ahora sí podemos pasar al kernel en CUDA C. Algunos nombres fueron cambiados debido al modo de escribir los programas, sin embargo la idea es la misma. Entre estos cambios notamos que blockDim fue cambiado a TAMANIO_BLOQUEy la introducción de tesela_idx que es en realidad blockDim*blockIdx + tx-1 con el que trabajamos en Python.
Este índice es usado puesto que para cada dato en A, necesitamos también su vecino a la izquierda (idx-1). Con tesela_idx cubrimos cada uno de estos. Ahora sólo falta el vecino derecho de la última celda. Este estará cubierto por otro pequeño programa condicional con un if else.
__shared__ float tesela_X[TAMANIO_BLOQUE+2] ;
int tx = threadIdx.x ;
int idx = blockIdx.x*TAMANIO_BLOQUE + tx ;
int tesela_idx = idx - 1 ;
if tx < TAMANIO_BLOQUE {
if ((tesela_idx >= 0) && (tesela_idx < Dim_Camino) ) {
tesela_X[tx] = X[tesela_idx] ;
} else {
tesela_X[tx] = 0.0f ;
}
__syncthreads() ;
}
if blockDim.x*(blockIdx.x+1) < Dim_Camino {
tesela_X[TAMANIO_BLOQUE+1] = X[blockDim.x*(blockIdx.x +1)] ;
} else {
tesela_X[TEMANIO_BLOQUE+1] = 0.0f ;
}
__syncthreads() ;
Una vez hecha la copia de los datos de X en tesela_X sólo falta reescribir X con los nuevos valores. Eso quedará de ustedes.
También es importante fijar el tamaño de los bloques TAMANIO_BLOQUE. Así que es hora de que el lector se ponga a trabajar y complete el kernel para luego graficar la evolución temporal de la distribución de probabilidad del caminante aleatorio. Supón en un primer tiempo que $p = q = \frac{1}{2}$.
Recomendamos graficar con la función imshow() de matplotlib.
A este método de resolver la ecuación maestra numéricamente se le llama enumeración exacta y es sumamente importante y utilizado para resolver ecuaciones diferenciales parciales.
[1b] Una vez que hayas obtenido las imágenes, cambia el valor de $p$ para ver como varía la distribución de probabilidad.
[2] Escribe la ecuación maestra del caminante aleatorio en dos dimensiones.
[3] Modifica tu código para obtener una seria de imágenes con las que puedas observar la evolución temporal de la distribución de probabilidad en 2 dimensiones.
Hint: En este caso tendrás que hacer una matriz tipo __shared__ y no un arreglo unidimensional. Te recomendamos revisar los notebooks sobre multiplicación de matrices para que recuerdes la manera de indexar.
Ahora tendremos cuatro indices:
int Fila = blockIdx.y*BLOCK_SIZE + ty ;
int Columna = blockIdx.x*BLOCK_SIZE + tx ;
int Fila_copia = Fila - 1 ; int Columna_copia = Columna - 1 ;
if( (Fila_copia >= 0) && (Fila_copia < DimY) && (Columna_copia >= 0) && (Columna_copia < DimX) ) {
ds_copiaPlano[ty][tx] = Plano[Fila_copia][Columna_copia];
} else {
ds_copiaPlano[ty][tx] = 0.0f ;
}
En caso de que te pierdas con los índices y copias, haz un código en Python semejante al de arriba que te sirva como guía.
Hasta ahora no nos hemos enfrentado con el problema de las fronteras, pero no podíamos escapar de él. Supongamos que son paredes reflejantes y no absorbentes, lo cual hará que el caminante "rebote" en las fronteras.
[4] Escribe la regla que tienen que seguir las probabilidades cuando un caminante llega a cualquiera de las cuatro fronteras.
[5] Implementa esta regla en tu código y observa qué pasa.
Supongamos ahora que los cambios en el espacio y en el tiempo del caminante aleatorio se dan por diferenciales $\delta x$ y $\delta t$. La ecuación (1) se vuelve entonces
$$ P(x, t+\delta t) = pP(x-\delta x, t) + qP(x+\delta x, t) $$Si ahora expandemos cada termino en series de Taylor (hasta 2º orden), llegamos a que:
$$ \frac{\partial P}{\partial t}(x, t) = (q-p)\frac{\delta x}{\delta t}\frac{\partial P}{\partial x}(x, t)+ \frac{\delta x^2}{2\delta t}\frac{\partial^2 P}{\partial x^2}(x, t) $$Si ahora volvemos de nuevo al caso $p = q = \frac{1}{2}$, y haciendo $D = \frac{\delta x^2}{2\delta t}$ obtenemos entonces la ya conocida ecuación de difusión.
$$\frac{\partial P}{\partial t}(x, t) = D\frac{\partial^2 P}{\partial x^2}(x, t)$$[6] Las soluciones analíticas de esta EDP son bien conocidas. Compáralas con tu solución numérica.
Así, vemos que el método de enumeración exacta para un caminante aleatorio provee un método numérico para resolver esta EDP de evolución. [El método se llama de diferencias finitas.]
In [ ]: